#### Packages ####library(targets)library(rnaturalearth)library(here)library(sf)library(gmRi)library(patchwork)library(gt)library(knitr)library(tidyverse)library(ggstream)library(ggforce)library(bcp)library(scales)library(corrplot)library(readxl)# Support functionssource(here("R/support/sizeSpectra_support.R"))# Set themetheme_set(theme_gmri())# Map polygonsus_poly <-ne_states("united states of america", returnclass ="sf")canada <-ne_states("canada", returnclass ="sf")# Function to process summaries for various factor combinationsget_group_summaries <-function(...){# Do some grouping to get totals group_totals <- nefsc_size_bins %>%group_by(...) %>%summarise(total_survey_catch =sum(numlen, na.rm = T),total_lw_bio =sum(sum_weight_kg, na.rm = T),total_strat_abund =sum(strat_total_abund_s, na.rm = T),total_strat_lw_bio =sum(strat_total_lwbio_s, na.rm = T), .groups ="drop") # length bins group_lengths <- nefsc_size_bins %>%group_by(..., length_bin) %>%summarise(lenbin_survey_catch =sum(numlen),lenbin_lw_bio =sum(sum_weight_kg),lenbin_strat_abund =sum(strat_total_abund_s),lenbin_strat_lw_bio =sum(strat_total_lwbio_s), .groups ="drop") %>%left_join(group_totals) %>%mutate(perc_total_catch = (lenbin_survey_catch - total_survey_catch) *100,perc_lw_bio = (lenbin_lw_bio - total_lw_bio) *100,perc_strat_abund = (lenbin_strat_abund - total_strat_abund) *100,perc_strat_lw_bio = (lenbin_strat_lw_bio - total_strat_lw_bio) *100)# weight bins group_weights <- nefsc_size_bins %>%group_by(..., weight_bin) %>%summarise(wtbin_survey_catch =sum(numlen),wtbin_lw_bio =sum(sum_weight_kg),wtbin_strat_abund =sum(strat_total_abund_s),wtbin_strat_lw_bio =sum(strat_total_lwbio_s),.groups ="drop") %>%left_join(group_totals) %>%mutate(perc_total_catch = (wtbin_survey_catch / total_survey_catch) *100,perc_lw_bio = (wtbin_lw_bio / total_lw_bio) *100,perc_strat_abund = (wtbin_strat_abund / total_strat_abund) *100,perc_strat_lw_bio = (wtbin_strat_lw_bio / total_strat_lw_bio) *100)return(list("length_bins"=drop_na(group_lengths),"weight_bins"=drop_na(group_weights)))}
Abstract
The Northwest Atlantic ocean bordering the United states and Scotian Shelf is one of the fastest warming locations on Earth. Research on the impacts of this rapid-warming has primarily focused on high-profile and/or upper trophic level species. Here we investigate an the ecosystem wide impacts by integrating community size structure changes using size spectrum analysis. While laboratory studies and ecological theory suggest that ectotherms raised at higher temperatures will reach smaller body-sizes at comparable development stages, it is unclear whether that relationship can be mitigated against across a diverse community through collective individual behaviors. In cases where communities fails to mitigate/adapt to elevated temperatures, we anticipate a steepening of the body-size spectrum slope relating to a reduction in larger sized individuals and an increase in smaller sized individuals within that community. Using groundfish abundance data from fisheries independent surveys we estimated size spectrum indices for four regions along the US NE continental shelf. At the regional scale we found that declines in the size spectrum slopes occurred the strongest in the Gulf of Maine and Georges Bank regions. Regions possessing the coldest temperatures and the largest populations of commercially targeted groundfish species. The largest declines in size spectrum slope occurred earlier in these areas than the more-recent temperature trends, suggesting other external forces may have driven the decline in larger-sized individuals within the communities. While the primary pressure of fisheries exploitation has declined over time, recovery in terms of larger-bodied individuals remains threatened by elevated temperatures.
Introduction
Temperature & Ecology
It is well understood that temperature plays a critical role on physiology through its impact on the chemical reactions that underpin biological life. A consequence of this, is that most life has evolved a thermal preference, around which the chemical reactions that support things like growth and means for self-preservation. Species that are unable to maintain their optimal thermal preferences internally (ectotherms), must be able to follow their thermal preferences through locomotion or adapt through changes in behavior. Otherwise; there will be physiological and/or metabolic costs in failing to do so.
In an era of anthropogenic climate change, there is an expectation that many species will be displaced from historic habitats. Research in marine environments has shown examples of species shifting to higher latitudes and deeper depths in the pursuit of more favorable conditions (Kleisner et al. 2017; Pinsky et al. 2013). Other research has suggested that the impacts of elevated temperatures may be manifesting not through geographic distribution change, but through physiological changes, changes in seasonal phenology, or hampered recovery efforts (Daan et al. 2005; Miller, O’Brien, and Fratantoni 2018; A. Pershing et al. 2018; University of South Carolina et al. 2021). Temperature has direct and indirect impacts on critical biological functions including the acquisition of biomass through feeding, the rates of biomass loss through metabolism, and the rates at which individuals mature and develop. This potential for temperature to impact the size structure of an ecosystem has implications for blue economy & natural resource systems we rely upon.
Size Spectrum Theory
Size is a defining characteristic of species that mediates many ecological interactions. Size impacts the mobility of an organism, its ability to evade predation, its ability to successfully prey on other organisms, and the metabolic costs associated with each of these behaviors. In the context of a strongly size-structured ecosystem, growth and maturity changes alter fitness and ultimately determine whether a species is successful in the given environment. Size structured environments are a fundamental organizational pattern that has been heavily researched. Ecological theory is rich with models relating how energy transfers from smaller prey species to larger predatory trophic levels, the allocation of energy for growth, and the trade offs of allocating energy towards reproduction at the expense of growth. One such ecological model that avoids the need for explicit articulation of each predator-prey interaction is the size spectrum model.
A “size spectrum” describes the distribution of biomass or abundance as a function of individuals’ mass or size on a log–log scale (Guiet, Poggiale, and Maury 2016). Size spectrum are described by two terms, the size spectrum slope & intercept. These two terms convey the baseline productivity, and how energy flows through an ecosystem in the form of biomass. The spectrum intercept value captures the richness or the productivity at the base of the community and is strongly connected to the prevailing environmental conditions (Boudreau and Dickie 1992). So much so, that spectrum shape is sometimes defined by its eutrophic or oligotrophic environmental conditions (Rossberg 2012).
Size spectra condense the complexities of predator prey networks and their interactions into a handful of size-based indices (SBI). In doing so a community size spectrum captures the emergent properties of a system, while needing only size and abundance information that is commonly available. For this reason it has become increasingly relied upon as an indicator of ecosystem health in the push for ecosystem based fisheries management. Changes in slope have been associated with fishing exploitation, primarily through the targeted removal of larger individuals (Bianchi et al. 2000; Shin et al. 2005). Numerical experiments have also tried to link changes in slope to environmental disturbances (Guiet, Poggiale, and Maury 2016). Biomass spectrum present predictable intercepts between ecosystems of similar productivity levels, but also of distinct temperatures (Guiet, Poggiale, and Maury 2016).
Because it controls chemical reactions, temperature controls metabolic rates which underpin maintenance, growth or reproduction (Clarke and Johnston, 1999; Kooijman, 2010) as well as the functional responses to food density (Rall et al., 2012). Guiet et al. (2016)… In addition to the impact of temperature on communities’ intercepts (heights), the impact of temperature on the speed of the energy flow within communities may affect other properties, such as their resilience to perturbations or the intensity of trophic cascades (Andersen and Pedersen, 2009).
Sea surface temperatures in the Gulf of Maine since 1982 have been warming at rates faster than 96% of the world’s oceans, with similar warming rates along the northwest Atlantic shelf (A. Pershing et al. 2018). A punctuated elevation in temperatures over the last decade are believed to be the result of a shift in Gulf Stream positioning. A Northward shift in the Gulf stream has increased the regional temperatures via an increased direct transport of warm gulf stream water into areas like the Gulf of Maine . The Gulf stream has also been producing a higher frequency of warm core rings, and has obstructed some of the cold-water Scotian Shelf current flow that would otherwise counter the influence of the Gulf Stream on the region’s temperatures (Gangopadhyay et al. 2019; University of South Carolina et al. 2021). The combination of these oceanographic changes has led to a warmer continental shelf habitat.
Species Trends in NE Shelf
The observed warming in the northwest Atlantic has been shown to be a major factor in the redistribution of marine species along the US east coast. Species have responded to these warming trends by adjusting the timing and locations of their seasonal migrations and shifting their geographic ranges (Nye et al. 2009; Staudinger et al. 2019). There is also evidence that warming has hampered fisheries recoveries through its physiological impacts on growth and metabolism. Species like Black Sea Bass, Atlantic shortfin squid, and Blue crab have been high-profile examples of species expanding their ranges to follow their thermal preferences. While species like the American lobster have shown declines at their southern range near Long Island Sound, with much doubt whether they will recover under the present temperature trends. The recent regime shift in the physical oceanography has also shown to be a catalyst for biological shifts as well (University of South Carolina et al. 2021; Perretti et al. 2017).
While these examples show that species can respond to changes in the physical environment around them through movement & behavior, research elsewhere suggests that physiological responses integrated across species will manifest as changes in community size structure.
Purpose
We estimated size spectrum relationships for the groundfish populations for each sub-region of the Northeast US continental shelf. In the case of the NW Atlantic: sustained increases in temperature should have a physiological impact on the community size structure.
This leads to our second hypothesis:
H2. Warming alters the community through the direct influence of temperature on metabolism, growth, and population productivity.
Methods
Groundfish Data
Code
# 1. Biological data used as inputwithr::with_dir(rprojroot::find_root('_targets.R'), tar_load(nefsc_1g_labelled)) # rename and formatnefsc_size_bins <- nefsc_1g_labelledn_species <-length(unique(nefsc_size_bins$comname))
Fishery Independent data on was collected as part of the NEFSC’s northeast trawl survey. This survey is conducted each year in the spring and in the fall, with sample locations determined following a stratified-random survey design with effort allocated in proportion to stratum area. Analyses using the trawl survey data took used data from both the Spring and Fall survey seasons covering all years from 1970 to 2019. Correction factors were applied to total species abundance and biomass to account for changes in vessels, gear, and doors when appropriate as part of the survey program. However, abundance and biomass by length is not corrected. To account for this, abundance at length for each species were adjusted to match the correction factors applied to total species abundance at each station, with allocation following the distributions of length caught at that station. Such that for each species: the sum of adjusted abundances for each length class is equal to the gross abundance as corrected for changes in vessels, gear, etc. To account for differences in sampling effort among survey strata, all corrected abundance-at-length data was area-stratified.
Code
# Get region polygons:# Get the paths to the shapefiles used to masktrawl_paths <- gmRi::get_timeseries_paths(region_group ="nmfs_trawl_regions", mac_os ="mojave")# Polygons for each regionall_poly <-read_sf(trawl_paths[["inuse_strata"]][["shape_path"]]) ggplot() +geom_sf(data = us_poly) +geom_sf(data = canada) +geom_sf(data = all_poly, aes(fill = survey_area), alpha =0.8) +coord_sf(xlim =c(-76.4, -64.4), ylim =c(35, 45.5), expand = F) +scale_fill_gmri() +theme_bw() +map_theme(legend.position ="bottom", legend.background =element_rect(fill ="white")) +guides(fill =guide_legend(title ="", nrow =2))
Data from the survey was grouped by strata to form geographically meaningful sub-regions: Gulf of Maine, Georges Bank, Southern New England, Mid-Atlantic Bight. For each region, we developed several time series indicators:
Community Composition metrics (abundance and biomass by functional group, with body-size contributions)
Mean size of the aggregate community and key functional groups
Slope and intercept of the size spectrum
Community Composition
Functional groups were assigned to each species based on life history and geography. Functional groups included were elasmobranch, pelagic, demersal, and groundfish species. Stratified annual abundance and biomass totals were calculated for each functional group and each region with labels for increasing body-size (biomass kg) groups.
Code
# make a column that contains all datanefsc_size_bins <- nefsc_size_bins %>%mutate(group_col ="All Data") # Do some grouping by year survey area and functional groupfgroup_area <-get_group_summaries(Year, survey_area, spec_class)
Body Size Trends
The annual stratified-abundance weighted body length and body weight within each region and for each functional group were also estimated using the numbers at length and the estimated biomass at length information. Data for body size trends was not truncated at 1g minimum body weight.
Code
# Load the average body size datawithr::with_dir(rprojroot::find_root('_targets.R'), tar_load(mean_sizes_ss_groups)) # Grouped on year and regionregional_size <- mean_sizes_ss_groups %>%filter(`group ID`=="single years * region") %>%mutate(Year =as.numeric(Year),factor(survey_area, levels =c("GoM", "GB", "SNE", "MAB")))
Size Spectrum Analysis
Community size spectra were estimated using abundance-at-length data from 68 species. These species were selected based on the availability of published weight-at-length relationships (Wigley et al. 2003) and represented 98.98% of the total biomass caught in the survey. Published length-weight relationships were used to convert abundance at length data into their corresponding biomass at length (kg). These values were then used to get totals for stratified weight-at-length, in complement to the corrected abundances-at-length data which had been area stratified. These area-stratified biomass at length totals were then used for fitting each regional biomass size spectra.
To fit the normalized biomass size spectra, stratified biomass at length data was binned into logartithmically equal spaced intervals (0.5 on a \(log_{10}\) scale), summing bodymass across all species within each body size bin. To normalize the spectra, the atratified abundances within each bin was then divided by the bin-width to account for the increasing bin-widths, a consequence of the log scale. Normalized size spectra were fit for each year and for each region independently, and for each year across all strata, using ordinary least squares (ols) regressions for stratified abundance (normalized) by body-size bins.
Code
withr::with_dir(rprojroot::find_root('_targets.R'), tar_load(size_spectrum_indices)) # Grab SS Groups we care aboutregion_indices <- size_spectrum_indices %>%filter(`group ID`=="single years * region") %>%mutate(yr =as.numeric(as.character(Year)),survey_area =factor(survey_area, levels =c("GoM", "GB", "SNE", "MAB")),sig_fit =ifelse(l10_sig_strat <0.05, "Significant", "Non-Significant"))
Temperature Data
Global Sea surface temperature data was obtained via NOAA’s optimally interpolated SST analysis (OISSTv2), providing daily temperature values at a 0.25° latitude x 0.25° longitude resolution (Reynolds et al. 2007). A daily climatology for every 0.25° pixel in the global data set was created using average daily temperatures spanning the period of 1982-2011. Daily anomalies were then computed as the difference between observed temperatures and the daily climatological average. OISSTv2 data used in these analyses were provided by the NOAA PSL, Boulder, Colorado, USA from their website at https://psl.noaa.gov.
Temperature data was regionally averaged to match the survey regions from the age-at-length data. SST anomalies were averaged by year for each region and over the entire sampling region to produce daily time series. These time series were then processed into annual timeseries of surface temperatures and anomalies. All region-averaging was done with area-weighting of the latitude/longitude grid cells to account for differences in cell-size in of the OISSTv2 data.
Code
# Load the regional temperatures from {targets}withr::with_dir(rprojroot::find_root('_targets.R'), tar_load(regional_oisst))
Spectra Drivers
The impact of external factors on the changes in size spectra was compared against several hypothesized driving forces. External factors proposed as potential drivers included sea surface temperature anomalies, state and federal fisheries landings from the Greater Atlantic Regional Fisheries Office (GARFO), the North Atlantic Oscillation (NAO), and zooplankton community indices from the continuous plankton recorder (CPR) dataset.
Code
#### 1. GARFO Landings ##### Load the GARFO landings data:res_path <- gmRi::cs_path("res")# Landings of finfish* sheet 5landings <-read_xlsx(path =str_c(res_path, "GARFO_landings/KMills_landings by area 1964-2021_JUN 2022.xlsx"), sheet =5) %>%rename_all(tolower)# Shapefiles for the fisheries stat zonesstat_zones <-read_sf(str_c(res_path, "Shapefiles/Statistical_Areas/Statistical_Areas_2010_withNames.shp"))# Make a list of zones to roughly match the survey areas:fish_zones <-list("Gulf of Maine"=c(511:515, 464, 465),"Georges Bank"=c(521, 522, 525, 561, 562),"Southern New England"=c(611, 612, 613, 616, 526, 537, 538, 539),"Mid-Atlantic Bight"=c(614:615, 621, 622, 625, 626, 631, 632))# Trim out what we don't need and label themstat_zones <- stat_zones %>%mutate(survey_area =case_when( Id %in% fish_zones$"Gulf of Maine"~"Gulf of Maine", Id %in% fish_zones$"Georges Bank"~"Georges Bank", Id %in% fish_zones$"Southern New England"~"Southern New England", Id %in% fish_zones$"Mid-Atlantic Bight"~"Mid-Atlantic Bight")) %>%filter(survey_area %in%c("Georges Bank", "Gulf of Maine", "Southern New England", "Mid-Atlantic Bight"))# # Map it out? - takes forever from Chesapeake bay# ggplot(st_simplify(stat_zones)) +# geom_sf(aes(fill = survey_area)) +# scale_fill_gmri() +# map_theme()# Add the labels into the landings data and remove what we don't need there:landings <- landings %>%mutate(survey_area =case_when(`stat area`%in% fish_zones$"Gulf of Maine"~"Gulf of Maine",`stat area`%in% fish_zones$"Georges Bank"~"Georges Bank",`stat area`%in% fish_zones$"Southern New England"~"Southern New England",`stat area`%in% fish_zones$"Mid-Atlantic Bight"~"Mid-Atlantic Bight")) %>%filter(survey_area %in%c("Georges Bank", "Gulf of Maine", "Southern New England", "Mid-Atlantic Bight")) %>%mutate(survey_area =factor(survey_area, levels =c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight" )))# Get Summarieslandings_summ <- landings %>%drop_na() %>%rename("weight_lb"=`landed lbs`,"live_lb"=`live lbs`) %>%drop_na() %>%group_by(year, survey_area) %>%summarise( across(.cols =c(value, weight_lb, live_lb), .fns =list(mean = mean, total = sum), .names ="{.fn}_{.col}"), .groups ="drop") # Scale the landings to create an index by arealandings_summ <- landings_summ %>%group_by(survey_area) %>%mutate(total_wt_z = base::scale(total_weight_lb)) %>%ungroup()#### 2. Climate Drivers ##### From ecodata:#### NAO & GSInao <- ecodata::nao %>%mutate(year = Time,Time =as.Date(str_c(Time, "-01-01")))gsi <- ecodata::gsi %>%mutate(Time =str_pad(as.character(Time), width =7, side ="right", pad ="0"),year =as.numeric(str_sub(Time, 1, 4)),month =str_sub(Time, -2, -1),Time =as.Date(str_c(year, month, "01", sep ="-")))# Join them upclim_drivers <-bind_rows(list(nao, gsi)) %>%filter(EPU =="All")#### 3. CPR Indices ##### From CPR Analysis:# Uses 7 focal taxa, 2 calanus stagescpr_indices <-read_csv("~/Documents/Repositories/continuous_plankton_recorder/results_data/cpr_focal_pca_timeseries_period_1961-2017.csv")# Do some Reshaping for the Plotcpr_indices <- cpr_indices %>%select(c(year, 2:3)) %>%rename(`Principal Component 1`="First Mode",`Principal Component 2`="Second Mode") %>%pivot_longer(names_to ="pca_mode", values_to ="pca_loading", cols =c(2:3))
Results
Abundance Distribution
Code
# panel plot the biomass by body sizefgroup_area$weight_bins %>%mutate(strat_abund_mill = wtbin_strat_abund /1e6) %>%ggplot(aes(Year, strat_abund_mill, fill = weight_bin)) +geom_col(position ="stack", color ="gray80", size =0.1) +facet_grid(survey_area~spec_class) +scale_fill_gmri() +scale_y_continuous(labels =comma_format()) +labs(y ="Stratified Total Abundance (millions)", title ="Abundance Allocation by Weight",x ="", fill ="Individual Weight (kg)") +theme(legend.position ="bottom",axis.text.x =element_text(angle =45, hjust =1))
Biomass Distribution
Overall biomass was highest in the two northern regions, the Gulf of Maine and Georges Bank. Roughly half of the biomass sampled in these regions can be attributed to groundfish/demersal species, with the second largest contributions coming from elasmobranchs. Groundfish biomass, larger individuals >2kg in particular, declined during the 70’s and 80’s in these regions, never truly recovering. Beginning in the 2000’s there were signs that groundfish abundances were increasing as evidenced by increasing numbers of smaller individuals, however in both regions this trend appears to have reversed by the mid 2010’s. Elasmobranch populations increased steadily throughout the survey time period across all regions, but for southern New England where they were more flat but with high/low periods. Larger elasmobranch were rare in all regions except for a period spanning the late 70’s through the early 90’s isolated to Georges Bank. Demersal species biomass was highest in the Gulf of Maine, dwarfing their contributions in other regions. Their biomass declined in the 70’s, was flat until the late 90’s, remaining relatively high until declining in the late 2010’s. Pelagic species biomass was low in all regions, and is unlikely to be representative of true biomass trends due to gear selectivity.
Code
# panel plot the biomass by body sizefgroup_area$weight_bins %>%mutate(strat_lwbio_mill = wtbin_strat_lw_bio /1e6) %>%ggplot(aes(Year, strat_lwbio_mill, fill = weight_bin)) +geom_col(position ="stack", color ="gray80", size =0.1) +facet_grid(survey_area~spec_class) +scale_fill_gmri() +scale_y_continuous(labels =comma_format()) +labs(y ="Stratified Total Biomass (million kg)", title ="Biomass Allocation by Weight",x ="", fill ="Individual Weight (kg)") +theme(legend.position ="bottom",axis.text.x =element_text(angle =45, hjust =1))
Regional Variation in Species Composition
There was a distinct difference between Northern and Southern regions in the way biomass was distributed among the different functional groups. The primary contributors to overall biomass in the southern regions (southern New England & mid-Atlantic bight) was the elasmobranch community. While the northern regions (Gulf of Maine & Georges Bank) each had similar quantities of elasmobranch biomasses, there was also a comparable contribution of groundfish and demersal species as well.
Size spectrum slopes were least steep in the two Northern regions.
Steepening of northern slopes
lack of bounce-back
all slopes similar steepness
Code
# Plot the Regional Slopesslope_timeline <-ggplot(region_indices, aes(yr, l10_slope_strat, group = survey_area)) +geom_line(linetype =1) +geom_point() +scale_color_gmri() +facet_wrap(~survey_area, ncol =1) +# geom_smooth(method = "lm", formula = y~x, alpha = 0.75) +labs(x ="Year", y ="Slope", title ="Normalized Biomass Spectra") # Plot the Regional Interceptsint_timeline <-ggplot(region_indices, aes(yr, l10_int_strat, group = survey_area)) +geom_line(linetype =1) +geom_point() +facet_wrap(~survey_area, ncol =1) +# geom_smooth(method = "lm", formula = y~x, alpha = 0.75) +labs(x ="Year", y ="Intercept") # Assembleslope_timeline | int_timeline
Size Spectra Drivers
Temperature Trends
Code
# Plot how each one exhibits a break in the temperature across the region# Get regional averagestemp_regimes <- regional_oisst %>%filter(yr >1981) %>%mutate(regime =ifelse(sst_anom >0, "hot", "cold"),survey_area =ifelse(survey_area =="all", "Full Region", survey_area),survey_area =factor(survey_area, levels =c("Full Region", "GoM", "GB", "SNE", "MAB")))
Overfishing Trends
Discussion
The data we rely on in this analysis was collected as part of a survey program which began out of concern that fisheries were already being over-harvested. Estimates by scientists at that time suggested that by the 1970’s total biomass of Georges Bank had been halved and elasmobranchs had begun to replace the over-exploited gadoid(Fogarty and Murawski 1998). The implication of such a large disturbance that pre-dates our time series is that the steepening of size spectrum slope we observed in this area and the adjacent Gulf of Maine are the tail-ends of a longer and more severe decline. While metrics of overall fishing pressure is hard to align exactly with trawl survey coverage, historic and anecdotal evidence show that groundfish fishing pressures are a fraction of their historic pressure was in the 1960’s and 1970’s. This begs the question of why larger adult numbers never began to recover in these regions. Looking at abundance and biomass information from the survey there was evidence of strong recruitment among smaller individuals < 1kg, but there has since not been any sizable population of fishes larger than 1kg outside of elasmobranchs. Work by (A. J. Pershing et al. 2015) suggests that part of the failure in recovery was due to an inability to account for temperature change in fisheries management. At this point in time the regional temperatures had just begun to reflect a regime shift, and could have been considered at that time an acute stressor. Since then the region has experienced nearly a decade of sustained above-average temperatures, and there are signs that the success seen in recruitment and survival of even the smaller size classes is declining.
#### CPR Community ##### CPR Plotcpr_indices %>%ggplot(aes(year, pca_loading, group = pca_mode)) +geom_hline(yintercept =0, color =gmri_cols("orange"), size =1, lty =2) +geom_line() +facet_wrap(~pca_mode, ncol =1) +scale_x_continuous(breaks =seq(1950, 2030, by =10)) +labs(x ="Year", y ="PCA Loading", title ="Gulf of Maine CPR PCA")
All Correlations Matrix
Code
# table to join for swapping shorthand for long-hand namesarea_df <-data.frame(area =c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "All"),survey_area =c("GoM", "GB", "SNE", "MAB", "Full Region"))# Common form for all:# year, area, var, value# Put climate drivers on annual scheduleclim_idx <- clim_drivers %>%group_by(year, area = EPU, var = Var) %>%summarise(value =mean(Value, na.rm = T),.groups ="drop")# landingslandings_idx <- landings_summ %>%group_by(year, area = survey_area) %>%summarise(value =mean(total_wt_z, na.rm = T),.groups ="drop") %>%mutate(var ="landings")# Zooplanktoncpr_idx <- cpr_indices %>%mutate(area ="Gulf of Maine",var =str_c("cpr_pca_", str_sub(pca_mode, -1, -1))) %>%group_by(year, area, var) %>%summarise(value =mean(pca_loading),.groups ="drop")# Temperaaturesst_idx <- temp_regimes %>%left_join(area_df) %>%select(year = yr, area, value = sst_anom) %>%mutate(var ="sst_anom")# Size Spectrum Slopes & Interceptsss_idx <- region_indices %>%select(year = Year, survey_area, spectra_slope = l10_slope_strat, spectra_int = l10_int_strat) %>%pivot_longer(cols =c(spectra_slope, spectra_int), names_to ="var", values_to ="value") %>%mutate(year =as.numeric(year)) %>%left_join(area_df) %>%select(-survey_area)# All driversall_drivers <-bind_rows(list( clim_idx, landings_idx, cpr_idx, sst_idx, ss_idx ))# Do some reshaping:driver_matrix <- all_drivers %>%mutate(id =str_replace_all(str_c(area, "_", var), " ", "_")) %>%select(-c(area, var)) %>%pivot_wider(names_from = id, values_from = value) %>%drop_na() %>%column_to_rownames(var ="year")# Correlation Plotrho <-cor(driver_matrix)corrplot(rho,method="ellipse",type="lower",order ="hclust",tl.col ="black",tl.srt =0,tl.cex =0.6,tl.offset =0.7,cl.cex =0.8,cl.offset =0.9,cl.ratio =0.1,title ="All Year Correlation Matrix")
gg_rho %>%mutate(xvar =str_replace(xvar, "_", " "),yvar =str_replace(yvar, "_", " ")) %>%ggplot(aes(yvar, xvar, fill = r2)) +geom_tile() +scale_fill_distiller(palette ="RdBu", limits =c(-1, 1), breaks =c(-1,0,1)) +theme(axis.text.x =element_text(angle =45, hjust =1, vjust =1),strip.text.y =element_text(angle =0),legend.position ="bottom") +labs(x ="With these Other Indices (y)", y ="Correlation of These Indices (x)",fill ="Correlation Coefficient",title ="Correlation Matrix of All Drivers & Response")
Code
# Tidy it up, enhance the labelsgg_rho_tidy <- gg_rho %>%mutate(response_var =ifelse(str_detect(xvar, "slope|int"), T, F),driver_var =ifelse(str_detect(yvar, "osc|pca|landings|index|anom"), T, F),area =case_when(str_detect(xvar, "Georges") ~"Georges Bank",str_detect(xvar, "Gulf") ~"Gulf of Maine",str_detect(xvar, "Southern") ~"Southern New England",str_detect(xvar, "Mid-Atlantic") ~"Mid-Atlantic Bight"),area =factor(area, levels =c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")),driver_type =ifelse(str_detect(yvar, "landings"), "Fishing", "Environmental"),x_short =case_when(str_detect(xvar, "slope") ~"Spectra Slope",str_detect(xvar, "int") ~"Spectra Intercept"))# Do some filtering and string formatting:only_driver_response <- gg_rho_tidy %>%filter( response_var, driver_var) %>%mutate(yvar =str_replace_all(yvar, "_", " "),yvar =factor(yvar, levels =c("Gulf of Maine cpr pca 1","Gulf of Maine cpr pca 2","All gulf stream index","All north atlantic oscillation","All sst anom","Gulf of Maine sst anom","Georges Bank sst anom", "Southern New England sst anom","Mid-Atlantic Bight sst anom","Gulf of Maine landings","Georges Bank landings","Southern New England landings","Mid-Atlantic Bight landings" ))) # Plot the convoluted thing:only_driver_response %>%ggplot(aes(yvar, x_short, fill = r2)) +geom_tile() +facet_grid(area~driver_type, scales ="free_x", space ="free") +scale_fill_distiller(palette ="RdBu", limits =c(-1, 1), breaks =c(-1,0,1)) +theme(axis.text.x =element_text(angle =45, hjust =1, vjust =1),strip.text.y =element_text(angle =0),legend.position ="bottom",legend.title =element_text(vjust =0.75)) +labs(x ="Driver", y ="Bodymass Spectra Coefficient",fill ="Correlation Coefficient",title ="Correlation Between Size Spectra Response and Drivers",subtitle ="Correlation for Concurrent Years")
References
Bianchi, G., H. Gislason, K. Graham, L. Hill, X. Jin, K. Koranteng, S. Manickchand-Heileman, et al. 2000. “Impact of Fishing on Size Composition and Diversity of Demersal Fish Communities.”ICES Journal of Marine Science 57 (3): 558–71. https://doi.org/10.1006/jmsc.2000.0727.
Boudreau, P. R., and L. M. Dickie. 1992. “Biomass Spectra of Aquatic Ecosystems in Relation to Fisheries Yield.”Canadian Journal of Fisheries and Aquatic Sciences 49 (8): 1528–38. https://doi.org/10.1139/f92-169.
Daan, Niels, Henrik Gislason, John G. Pope, and Jake C. Rice. 2005. “Changes in the North Sea Fish Community: Evidence of Indirect Effects of Fishing?”ICES Journal of Marine Science 62 (2): 177–88. https://doi.org/10.1016/j.icesjms.2004.08.020.
Fogarty, Michael J., and Steven A. Murawski. 1998. “Large-Scale Disturbance and the Structure of Marine Systems: Fishery Impacts on Georges Bank.”Ecological Applications 8 (sp1): S6–22. https://doi.org/10.1890/1051-0761(1998)8[S6:LDATSO]2.0.CO;2.
Gangopadhyay, Avijit, Glen Gawarkiewicz, E. Nishchitha S. Silva, M. Monim, and Jenifer Clark. 2019. “An Observed Regime Shift in the Formation of Warm Core Rings from the Gulf Stream.”Scientific Reports 9 (1): 12319. https://doi.org/10.1038/s41598-019-48661-9.
Guiet, Jérôme, Jean-Christophe Poggiale, and Olivier Maury. 2016. “Modelling the Community Size-Spectrum: Recent Developments and New Directions.”Ecological Modelling 337 (October): 4–14. https://doi.org/10.1016/j.ecolmodel.2016.05.015.
Kleisner, Kristin M., Michael J. Fogarty, Sally McGee, Jonathan A. Hare, Skye Moret, Charles T. Perretti, and Vincent S. Saba. 2017. “Marine Species Distribution Shifts on the U.S. Northeast Continental Shelf Under Continued Ocean Warming.”Progress in Oceanography 153 (April): 24–36. https://doi.org/10.1016/j.pocean.2017.04.001.
Miller, Timothy J., Loretta O’Brien, and Paula S. Fratantoni. 2018. “Temporal and Environmental Variation in Growth and Maturity and Effects on Management Reference Points of Georges Bank Atlantic Cod.”Canadian Journal of Fisheries and Aquatic Sciences 75 (12): 2159–71. https://doi.org/10.1139/cjfas-2017-0124.
Nye, Ja, Js Link, Ja Hare, and Wj Overholtz. 2009. “Changing Spatial Distribution of Fish Stocks in Relation to Climate and Population Size on the Northeast United States Continental Shelf.”Marine Ecology Progress Series 393 (October): 111–29. https://doi.org/10.3354/meps08220.
Perretti, Ct, Mj Fogarty, Kd Friedland, Ja Hare, Sm Lucey, Rs McBride, Tj Miller, et al. 2017. “Regime Shifts in Fish Recruitment on the Northeast US Continental Shelf.”Marine Ecology Progress Series 574 (July): 1–11. https://doi.org/10.3354/meps12183.
Pershing, Andrew J., Michael A. Alexander, Christina M. Hernandez, Lisa A. Kerr, Arnault Le Bris, Katherine E. Mills, Janet A. Nye, et al. 2015. “Slow Adaptation in the Face of Rapid Warming Leads to Collapse of the Gulf of Maine Cod Fishery.”Science 350 (6262): 809–12. https://doi.org/10.1126/science.aac9819.
Pershing, Andrew, Katherine Mills, Alexa Dayton, Bradley Franklin, and Brian Kennedy. 2018. “Evidence for Adaptation from the 2016 Marine Heatwave in the Northwest Atlantic Ocean.”Oceanography 31 (2). https://doi.org/10.5670/oceanog.2018.213.
Pinsky, Malin L., Boris Worm, Michael J. Fogarty, Jorge L. Sarmiento, and Simon A. Levin. 2013. “Marine Taxa Track Local Climate Velocities.”Science 341 (6151): 1239–42. https://doi.org/10.1126/science.1239352.
Rossberg, Axel G. 2012. “6 - A Complete Analytic Theory for Structure and Dynamics of Populations and Communities Spanning Wide Ranges in Body Size.” In, edited by Ute Jacob and Guy Woodward, 46:427–521. Global Change in Multispecies Systems Part 1. Academic Press. https://doi.org/10.1016/B978-0-12-396992-7.00008-3.
Shin, Yunne-Jai, Marie-Joëlle Rochet, Simon Jennings, John G. Field, and Henrik Gislason. 2005. “Using Size-Based Indicators to Evaluate the Ecosystem Effects of Fishing.”ICES Journal of Marine Science 62 (3): 384–96. https://doi.org/10.1016/j.icesjms.2005.01.004.
Staudinger, Michelle D., Katherine E. Mills, Karen Stamieszkin, Nicholas R. Record, Christine A. Hudak, Andrew Allyn, Antony Diamond, et al. 2019. “It’s about Time: A Synthesis of Changing Phenology in the Gulf of Maine Ecosystem.”Fisheries Oceanography 28 (5): 532–66. https://doi.org/10.1111/fog.12429.
University of South Carolina, Erin Meyer-Gutbrod, Charles Greene, Kimberley Davies, and David Johns. 2021. “Ocean Regime Shift Is Driving Collapse of the North Atlantic Right Whale Population.”Oceanography 34 (3): 22–31. https://doi.org/10.5670/oceanog.2021.308.
Source Code
---title: "Temperature Impacts on Size Structure of Northeast US Groundfish Community Extend Issues from Fisheries Over-Exploitation"author: name: "Adam Kemberling" url: https://github.com/adamkemberling affiliation: Gulf of Maine Research Institutedescription: | Size Spectrum Analysis of the Northeast US Groundfish Communitydate: "`r Sys.Date()`"format: html: self-contained: true code-fold: true code-tools: true df-print: kable toc: true toc-depth: 2editor: sourceexecute: echo: true warning: false message: false fig.height: 6 fig.width: 6 fig.align: "center" comment: ""bibliography: references.bib---```{r}#| label: load packages and functions#### Packages ####library(targets)library(rnaturalearth)library(here)library(sf)library(gmRi)library(patchwork)library(gt)library(knitr)library(tidyverse)library(ggstream)library(ggforce)library(bcp)library(scales)library(corrplot)library(readxl)# Support functionssource(here("R/support/sizeSpectra_support.R"))# Set themetheme_set(theme_gmri())# Map polygonsus_poly <-ne_states("united states of america", returnclass ="sf")canada <-ne_states("canada", returnclass ="sf")# Function to process summaries for various factor combinationsget_group_summaries <-function(...){# Do some grouping to get totals group_totals <- nefsc_size_bins %>%group_by(...) %>%summarise(total_survey_catch =sum(numlen, na.rm = T),total_lw_bio =sum(sum_weight_kg, na.rm = T),total_strat_abund =sum(strat_total_abund_s, na.rm = T),total_strat_lw_bio =sum(strat_total_lwbio_s, na.rm = T), .groups ="drop") # length bins group_lengths <- nefsc_size_bins %>%group_by(..., length_bin) %>%summarise(lenbin_survey_catch =sum(numlen),lenbin_lw_bio =sum(sum_weight_kg),lenbin_strat_abund =sum(strat_total_abund_s),lenbin_strat_lw_bio =sum(strat_total_lwbio_s), .groups ="drop") %>%left_join(group_totals) %>%mutate(perc_total_catch = (lenbin_survey_catch - total_survey_catch) *100,perc_lw_bio = (lenbin_lw_bio - total_lw_bio) *100,perc_strat_abund = (lenbin_strat_abund - total_strat_abund) *100,perc_strat_lw_bio = (lenbin_strat_lw_bio - total_strat_lw_bio) *100)# weight bins group_weights <- nefsc_size_bins %>%group_by(..., weight_bin) %>%summarise(wtbin_survey_catch =sum(numlen),wtbin_lw_bio =sum(sum_weight_kg),wtbin_strat_abund =sum(strat_total_abund_s),wtbin_strat_lw_bio =sum(strat_total_lwbio_s),.groups ="drop") %>%left_join(group_totals) %>%mutate(perc_total_catch = (wtbin_survey_catch / total_survey_catch) *100,perc_lw_bio = (wtbin_lw_bio / total_lw_bio) *100,perc_strat_abund = (wtbin_strat_abund / total_strat_abund) *100,perc_strat_lw_bio = (wtbin_strat_lw_bio / total_strat_lw_bio) *100)return(list("length_bins"=drop_na(group_lengths),"weight_bins"=drop_na(group_weights)))}```# AbstractThe Northwest Atlantic ocean bordering the United states and Scotian Shelf is one of the fastest warming locations on Earth. Research on the impacts of this rapid-warming has primarily focused on high-profile and/or upper trophic level species. Here we investigate an the ecosystem wide impacts by integrating community size structure changes using size spectrum analysis. While laboratory studies and ecological theory suggest that ectotherms raised at higher temperatures will reach smaller body-sizes at comparable development stages, it is unclear whether that relationship can be mitigated against across a diverse community through collective individual behaviors. In cases where communities fails to mitigate/adapt to elevated temperatures, we anticipate a steepening of the body-size spectrum slope relating to a reduction in larger sized individuals and an increase in smaller sized individuals within that community. Using groundfish abundance data from fisheries independent surveys we estimated size spectrum indices for four regions along the US NE continental shelf. At the regional scale we found that declines in the size spectrum slopes occurred the strongest in the Gulf of Maine and Georges Bank regions. Regions possessing the coldest temperatures and the largest populations of commercially targeted groundfish species. The largest declines in size spectrum slope occurred earlier in these areas than the more-recent temperature trends, suggesting other external forces may have driven the decline in larger-sized individuals within the communities. While the primary pressure of fisheries exploitation has declined over time, recovery in terms of larger-bodied individuals remains threatened by elevated temperatures.# Introduction## Temperature & EcologyIt is well understood that temperature plays a critical role on physiology through its impact on the chemical reactions that underpin biological life. A consequence of this, is that most life has evolved a thermal preference, around which the chemical reactions that support things like growth and means for self-preservation. Species that are unable to maintain their optimal thermal preferences internally (ectotherms), must be able to follow their thermal preferences through locomotion or adapt through changes in behavior. Otherwise; there will be physiological and/or metabolic costs in failing to do so.In an era of anthropogenic climate change, there is an expectation that many species will be displaced from historic habitats. Research in marine environments has shown examples of species shifting to higher latitudes and deeper depths in the pursuit of more favorable conditions [@kleisner2017; @pinsky2013]. Other research has suggested that the impacts of elevated temperatures may be manifesting not through geographic distribution change, but through physiological changes, changes in seasonal phenology, or hampered recovery efforts [@daan2005; @miller2018; @pershing2018; @meyer-gutbrod2021]. Temperature has direct and indirect impacts on critical biological functions including the acquisition of biomass through feeding, the rates of biomass loss through metabolism, and the rates at which individuals mature and develop. This potential for temperature to impact the size structure of an ecosystem has implications for blue economy & natural resource systems we rely upon.## Size Spectrum TheorySize is a defining characteristic of species that mediates many ecological interactions. Size impacts the mobility of an organism, its ability to evade predation, its ability to successfully prey on other organisms, and the metabolic costs associated with each of these behaviors. In the context of a strongly size-structured ecosystem, growth and maturity changes alter fitness and ultimately determine whether a species is successful in the given environment. Size structured environments are a fundamental organizational pattern that has been heavily researched. Ecological theory is rich with models relating how energy transfers from smaller prey species to larger predatory trophic levels, the allocation of energy for growth, and the trade offs of allocating energy towards reproduction at the expense of growth. One such ecological model that avoids the need for explicit articulation of each predator-prey interaction is the size spectrum model.A "size spectrum" describes the distribution of biomass or abundance as a function of individuals' mass or size on a log--log scale [@guiet2016]. Size spectrum are described by two terms, the size spectrum slope & intercept. These two terms convey the baseline productivity, and how energy flows through an ecosystem in the form of biomass. The spectrum intercept value captures the richness or the productivity at the base of the community and is strongly connected to the prevailing environmental conditions [@boudreau1992]. So much so, that spectrum shape is sometimes defined by its eutrophic or oligotrophic environmental conditions [@rossberg2012].Size spectra condense the complexities of predator prey networks and their interactions into a handful of size-based indices (SBI). In doing so a community size spectrum captures the emergent properties of a system, while needing only size and abundance information that is commonly available. For this reason it has become increasingly relied upon as an indicator of ecosystem health in the push for ecosystem based fisheries management. Changes in slope have been associated with fishing exploitation, primarily through the targeted removal of larger individuals [@bianchi2000; @shin2005]. Numerical experiments have also tried to link changes in slope to environmental disturbances [@guiet2016]. Biomass spectrum present predictable intercepts between ecosystems of similar productivity levels, but also of distinct temperatures [@guiet2016].This is a direct quote from [@guiet2016], but nails the connection back to temp expectations:> Because it controls chemical reactions, temperature controls metabolic rates which underpin maintenance, growth or reproduction (Clarke and Johnston, 1999; Kooijman, 2010) as well as the functional responses to food density (Rall et al., 2012). Guiet et al. (2016)... In addition to the impact of temperature on communities' intercepts (heights), the impact of temperature on the speed of the energy flow within communities may affect other properties, such as their resilience to perturbations or the intensity of trophic cascades (Andersen and Pedersen, 2009).### Temperature of the Gulf of Maine & NE Shelf- Marine species distribution shifts in NE Atlantic shelf [@kleisner2017]Sea surface temperatures in the Gulf of Maine since 1982 have been warming at rates faster than 96% of the world's oceans, with similar warming rates along the northwest Atlantic shelf [@pershing2018]. A punctuated elevation in temperatures over the last decade are believed to be the result of a shift in Gulf Stream positioning. A Northward shift in the Gulf stream has increased the regional temperatures via an increased direct transport of warm gulf stream water into areas like the Gulf of Maine . The Gulf stream has also been producing a higher frequency of warm core rings, and has obstructed some of the cold-water Scotian Shelf current flow that would otherwise counter the influence of the Gulf Stream on the region's temperatures [@gangopadhyay2019; @meyer-gutbrod2021]. The combination of these oceanographic changes has led to a warmer continental shelf habitat.### Species Trends in NE ShelfThe observed warming in the northwest Atlantic has been shown to be a major factor in the redistribution of marine species along the US east coast. Species have responded to these warming trends by adjusting the timing and locations of their seasonal migrations and shifting their geographic ranges [@nye2009; @staudinger2019]. There is also evidence that warming has hampered fisheries recoveries through its physiological impacts on growth and metabolism. Species like Black Sea Bass, Atlantic shortfin squid, and Blue crab have been high-profile examples of species expanding their ranges to follow their thermal preferences. While species like the American lobster have shown declines at their southern range near Long Island Sound, with much doubt whether they will recover under the present temperature trends. The recent regime shift in the physical oceanography has also shown to be a catalyst for biological shifts as well [@meyer-gutbrod2021; @perretti2017].While these examples show that species can respond to changes in the physical environment around them through movement & behavior, research elsewhere suggests that physiological responses integrated across species will manifest as changes in community size structure.### PurposeWe estimated size spectrum relationships for the groundfish populations for each sub-region of the Northeast US continental shelf. In the case of the NW Atlantic: sustained increases in temperature should have a physiological impact on the community size structure.This leads to our second hypothesis:> #### H2. Warming alters the community through the direct influence of temperature on metabolism, growth, and population productivity.# Methods## Groundfish Data```{r}#| label: load survdat_data# 1. Biological data used as inputwithr::with_dir(rprojroot::find_root('_targets.R'), tar_load(nefsc_1g_labelled)) # rename and formatnefsc_size_bins <- nefsc_1g_labelledn_species <-length(unique(nefsc_size_bins$comname))```Fishery Independent data on was collected as part of the NEFSC's northeast trawl survey. This survey is conducted each year in the spring and in the fall, with sample locations determined following a stratified-random survey design with effort allocated in proportion to stratum area. Analyses using the trawl survey data took used data from both the Spring and Fall survey seasons covering all years from 1970 to 2019. Correction factors were applied to total species abundance and biomass to account for changes in vessels, gear, and doors when appropriate as part of the survey program. However, abundance and biomass by length is not corrected. To account for this, abundance at length for each species were adjusted to match the correction factors applied to total species abundance at each station, with allocation following the distributions of length caught at that station. Such that for each species: the sum of adjusted abundances for each length class is equal to the gross abundance as corrected for changes in vessels, gear, etc. To account for differences in sampling effort among survey strata, all corrected abundance-at-length data was area-stratified.```{r}#| label: make region map# Get region polygons:# Get the paths to the shapefiles used to masktrawl_paths <- gmRi::get_timeseries_paths(region_group ="nmfs_trawl_regions", mac_os ="mojave")# Polygons for each regionall_poly <-read_sf(trawl_paths[["inuse_strata"]][["shape_path"]]) ggplot() +geom_sf(data = us_poly) +geom_sf(data = canada) +geom_sf(data = all_poly, aes(fill = survey_area), alpha =0.8) +coord_sf(xlim =c(-76.4, -64.4), ylim =c(35, 45.5), expand = F) +scale_fill_gmri() +theme_bw() +map_theme(legend.position ="bottom", legend.background =element_rect(fill ="white")) +guides(fill =guide_legend(title ="", nrow =2))```Data from the survey was grouped by strata to form geographically meaningful sub-regions: Gulf of Maine, Georges Bank, Southern New England, Mid-Atlantic Bight. For each region, we developed several time series indicators:1. Community Composition metrics (abundance and biomass by functional group, with body-size contributions)2. Mean size of the aggregate community and key functional groups3. Slope and intercept of the size spectrum\### Community CompositionFunctional groups were assigned to each species based on life history and geography. Functional groups included were elasmobranch, pelagic, demersal, and groundfish species. Stratified annual abundance and biomass totals were calculated for each functional group and each region with labels for increasing body-size (biomass kg) groups.```{r}#| label: group biomass and abundance data# make a column that contains all datanefsc_size_bins <- nefsc_size_bins %>%mutate(group_col ="All Data") # Do some grouping by year survey area and functional groupfgroup_area <-get_group_summaries(Year, survey_area, spec_class)```### Body Size TrendsThe annual stratified-abundance weighted body length and body weight within each region and for each functional group were also estimated using the numbers at length and the estimated biomass at length information. Data for body size trends was not truncated at 1g minimum body weight.```{r}#| label: body size data# Load the average body size datawithr::with_dir(rprojroot::find_root('_targets.R'), tar_load(mean_sizes_ss_groups)) # Grouped on year and regionregional_size <- mean_sizes_ss_groups %>%filter(`group ID`=="single years * region") %>%mutate(Year =as.numeric(Year),factor(survey_area, levels =c("GoM", "GB", "SNE", "MAB")))```### Size Spectrum AnalysisCommunity size spectra were estimated using abundance-at-length data from `r n_species` species. These species were selected based on the availability of published weight-at-length relationships (Wigley et al. 2003) and represented 98.98% of the total biomass caught in the survey. Published length-weight relationships were used to convert abundance at length data into their corresponding biomass at length (kg). These values were then used to get totals for stratified weight-at-length, in complement to the corrected abundances-at-length data which had been area stratified. These area-stratified biomass at length totals were then used for fitting each regional biomass size spectra.To fit the normalized biomass size spectra, stratified biomass at length data was binned into logartithmically equal spaced intervals (0.5 on a $log_{10}$ scale), summing bodymass across all species within each body size bin. To normalize the spectra, the atratified abundances within each bin was then divided by the bin-width to account for the increasing bin-widths, a consequence of the log scale. Normalized size spectra were fit for each year and for each region independently, and for each year across all strata, using ordinary least squares (ols) regressions for stratified abundance (normalized) by body-size bins.```{r}#| label: load the size spectrum indiceswithr::with_dir(rprojroot::find_root('_targets.R'), tar_load(size_spectrum_indices)) # Grab SS Groups we care aboutregion_indices <- size_spectrum_indices %>%filter(`group ID`=="single years * region") %>%mutate(yr =as.numeric(as.character(Year)),survey_area =factor(survey_area, levels =c("GoM", "GB", "SNE", "MAB")),sig_fit =ifelse(l10_sig_strat <0.05, "Significant", "Non-Significant"))```## Temperature DataGlobal Sea surface temperature data was obtained via NOAA's optimally interpolated SST analysis (OISSTv2), providing daily temperature values at a 0.25° latitude x 0.25° longitude resolution (Reynolds et al. 2007). A daily climatology for every 0.25° pixel in the global data set was created using average daily temperatures spanning the period of 1982-2011. Daily anomalies were then computed as the difference between observed temperatures and the daily climatological average. OISSTv2 data used in these analyses were provided by the NOAA PSL, Boulder, Colorado, USA from their website at https://psl.noaa.gov.Temperature data was regionally averaged to match the survey regions from the age-at-length data. SST anomalies were averaged by year for each region and over the entire sampling region to produce daily time series. These time series were then processed into annual timeseries of surface temperatures and anomalies. All region-averaging was done with area-weighting of the latitude/longitude grid cells to account for differences in cell-size in of the OISSTv2 data.```{r}#| label: load temperature data# Load the regional temperatures from {targets}withr::with_dir(rprojroot::find_root('_targets.R'), tar_load(regional_oisst))```## Spectra DriversThe impact of external factors on the changes in size spectra was compared against several hypothesized driving forces. External factors proposed as potential drivers included sea surface temperature anomalies, state and federal fisheries landings from the Greater Atlantic Regional Fisheries Office (GARFO), the North Atlantic Oscillation (NAO), and zooplankton community indices from the continuous plankton recorder (CPR) dataset.```{r}#| label: fishing-effort#| eval: true#| echo: true#### 1. GARFO Landings ##### Load the GARFO landings data:res_path <- gmRi::cs_path("res")# Landings of finfish* sheet 5landings <-read_xlsx(path =str_c(res_path, "GARFO_landings/KMills_landings by area 1964-2021_JUN 2022.xlsx"), sheet =5) %>%rename_all(tolower)# Shapefiles for the fisheries stat zonesstat_zones <-read_sf(str_c(res_path, "Shapefiles/Statistical_Areas/Statistical_Areas_2010_withNames.shp"))# Make a list of zones to roughly match the survey areas:fish_zones <-list("Gulf of Maine"=c(511:515, 464, 465),"Georges Bank"=c(521, 522, 525, 561, 562),"Southern New England"=c(611, 612, 613, 616, 526, 537, 538, 539),"Mid-Atlantic Bight"=c(614:615, 621, 622, 625, 626, 631, 632))# Trim out what we don't need and label themstat_zones <- stat_zones %>%mutate(survey_area =case_when( Id %in% fish_zones$"Gulf of Maine"~"Gulf of Maine", Id %in% fish_zones$"Georges Bank"~"Georges Bank", Id %in% fish_zones$"Southern New England"~"Southern New England", Id %in% fish_zones$"Mid-Atlantic Bight"~"Mid-Atlantic Bight")) %>%filter(survey_area %in%c("Georges Bank", "Gulf of Maine", "Southern New England", "Mid-Atlantic Bight"))# # Map it out? - takes forever from Chesapeake bay# ggplot(st_simplify(stat_zones)) +# geom_sf(aes(fill = survey_area)) +# scale_fill_gmri() +# map_theme()# Add the labels into the landings data and remove what we don't need there:landings <- landings %>%mutate(survey_area =case_when(`stat area`%in% fish_zones$"Gulf of Maine"~"Gulf of Maine",`stat area`%in% fish_zones$"Georges Bank"~"Georges Bank",`stat area`%in% fish_zones$"Southern New England"~"Southern New England",`stat area`%in% fish_zones$"Mid-Atlantic Bight"~"Mid-Atlantic Bight")) %>%filter(survey_area %in%c("Georges Bank", "Gulf of Maine", "Southern New England", "Mid-Atlantic Bight")) %>%mutate(survey_area =factor(survey_area, levels =c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight" )))# Get Summarieslandings_summ <- landings %>%drop_na() %>%rename("weight_lb"=`landed lbs`,"live_lb"=`live lbs`) %>%drop_na() %>%group_by(year, survey_area) %>%summarise( across(.cols =c(value, weight_lb, live_lb), .fns =list(mean = mean, total = sum), .names ="{.fn}_{.col}"), .groups ="drop") # Scale the landings to create an index by arealandings_summ <- landings_summ %>%group_by(survey_area) %>%mutate(total_wt_z = base::scale(total_weight_lb)) %>%ungroup()#### 2. Climate Drivers ##### From ecodata:#### NAO & GSInao <- ecodata::nao %>%mutate(year = Time,Time =as.Date(str_c(Time, "-01-01")))gsi <- ecodata::gsi %>%mutate(Time =str_pad(as.character(Time), width =7, side ="right", pad ="0"),year =as.numeric(str_sub(Time, 1, 4)),month =str_sub(Time, -2, -1),Time =as.Date(str_c(year, month, "01", sep ="-")))# Join them upclim_drivers <-bind_rows(list(nao, gsi)) %>%filter(EPU =="All")#### 3. CPR Indices ##### From CPR Analysis:# Uses 7 focal taxa, 2 calanus stagescpr_indices <-read_csv("~/Documents/Repositories/continuous_plankton_recorder/results_data/cpr_focal_pca_timeseries_period_1961-2017.csv")# Do some Reshaping for the Plotcpr_indices <- cpr_indices %>%select(c(year, 2:3)) %>%rename(`Principal Component 1`="First Mode",`Principal Component 2`="Second Mode") %>%pivot_longer(names_to ="pca_mode", values_to ="pca_loading", cols =c(2:3)) ```# Results### Abundance Distribution```{r}#| label: abundance distributions#| fig.width: 8#| fig-height: 8# panel plot the biomass by body sizefgroup_area$weight_bins %>%mutate(strat_abund_mill = wtbin_strat_abund /1e6) %>%ggplot(aes(Year, strat_abund_mill, fill = weight_bin)) +geom_col(position ="stack", color ="gray80", size =0.1) +facet_grid(survey_area~spec_class) +scale_fill_gmri() +scale_y_continuous(labels =comma_format()) +labs(y ="Stratified Total Abundance (millions)", title ="Abundance Allocation by Weight",x ="", fill ="Individual Weight (kg)") +theme(legend.position ="bottom",axis.text.x =element_text(angle =45, hjust =1))```### Biomass DistributionOverall biomass was highest in the two northern regions, the Gulf of Maine and Georges Bank. Roughly half of the biomass sampled in these regions can be attributed to groundfish/demersal species, with the second largest contributions coming from elasmobranchs. Groundfish biomass, larger individuals \>2kg in particular, declined during the 70's and 80's in these regions, never truly recovering. Beginning in the 2000's there were signs that groundfish abundances were increasing as evidenced by increasing numbers of smaller individuals, however in both regions this trend appears to have reversed by the mid 2010's. Elasmobranch populations increased steadily throughout the survey time period across all regions, but for southern New England where they were more flat but with high/low periods. Larger elasmobranch were rare in all regions except for a period spanning the late 70's through the early 90's isolated to Georges Bank. Demersal species biomass was highest in the Gulf of Maine, dwarfing their contributions in other regions. Their biomass declined in the 70's, was flat until the late 90's, remaining relatively high until declining in the late 2010's. Pelagic species biomass was low in all regions, and is unlikely to be representative of true biomass trends due to gear selectivity.```{r}#| label: bodymass distributions#| fig.width: 8#| fig-height: 8# panel plot the biomass by body sizefgroup_area$weight_bins %>%mutate(strat_lwbio_mill = wtbin_strat_lw_bio /1e6) %>%ggplot(aes(Year, strat_lwbio_mill, fill = weight_bin)) +geom_col(position ="stack", color ="gray80", size =0.1) +facet_grid(survey_area~spec_class) +scale_fill_gmri() +scale_y_continuous(labels =comma_format()) +labs(y ="Stratified Total Biomass (million kg)", title ="Biomass Allocation by Weight",x ="", fill ="Individual Weight (kg)") +theme(legend.position ="bottom",axis.text.x =element_text(angle =45, hjust =1))```### Regional Variation in Species CompositionThere was a distinct difference between Northern and Southern regions in the way biomass was distributed among the different functional groups. The primary contributors to overall biomass in the southern regions (southern New England & mid-Atlantic bight) was the elasmobranch community. While the northern regions (Gulf of Maine & Georges Bank) each had similar quantities of elasmobranch biomasses, there was also a comparable contribution of groundfish and demersal species as well.### Body Size Trends```{r}#| label: average body size trends#| fig.height: 8# Length plotavg_len_p <-ggplot(regional_size, aes(Year, mean_len_cm, group = survey_area)) +geom_line() +geom_point() +facet_wrap(~survey_area, ncol =1) +labs(title ="Average Length",y ="Length (cm)")# Weight plotavg_wt_p <-ggplot(regional_size, aes(Year, mean_wt_kg, group = survey_area)) +geom_line() +geom_point() +facet_wrap(~survey_area, ncol =1)+labs(title ="Average Weight",y ="Weight (kg)")# Plot side by sideavg_len_p | avg_wt_p```### Regional Size SpectraSize spectrum slopes were least steep in the two Northern regions.- Steepening of northern slopes- lack of bounce-back- all slopes similar steepness```{r}#| label: size spectra results#| fig.height: 8# Plot the Regional Slopesslope_timeline <-ggplot(region_indices, aes(yr, l10_slope_strat, group = survey_area)) +geom_line(linetype =1) +geom_point() +scale_color_gmri() +facet_wrap(~survey_area, ncol =1) +# geom_smooth(method = "lm", formula = y~x, alpha = 0.75) +labs(x ="Year", y ="Slope", title ="Normalized Biomass Spectra") # Plot the Regional Interceptsint_timeline <-ggplot(region_indices, aes(yr, l10_int_strat, group = survey_area)) +geom_line(linetype =1) +geom_point() +facet_wrap(~survey_area, ncol =1) +# geom_smooth(method = "lm", formula = y~x, alpha = 0.75) +labs(x ="Year", y ="Intercept") # Assembleslope_timeline | int_timeline ```## Size Spectra Drivers### Temperature Trends```{r}#| label: Temperature regime shift#| fig.height: 8# Plot how each one exhibits a break in the temperature across the region# Get regional averagestemp_regimes <- regional_oisst %>%filter(yr >1981) %>%mutate(regime =ifelse(sst_anom >0, "hot", "cold"),survey_area =ifelse(survey_area =="all", "Full Region", survey_area),survey_area =factor(survey_area, levels =c("Full Region", "GoM", "GB", "SNE", "MAB"))) ```### Overfishing Trends# DiscussionThe data we rely on in this analysis was collected as part of a survey program which began out of concern that fisheries were already being over-harvested. Estimates by scientists at that time suggested that by the 1970's total biomass of Georges Bank had been halved and elasmobranchs had begun to replace the over-exploited gadoid[@fogarty1998]. The implication of such a large disturbance that pre-dates our time series is that the steepening of size spectrum slope we observed in this area and the adjacent Gulf of Maine are the tail-ends of a longer and more severe decline. While metrics of overall fishing pressure is hard to align exactly with trawl survey coverage, historic and anecdotal evidence show that groundfish fishing pressures are a fraction of their historic pressure was in the 1960's and 1970's. This begs the question of why larger adult numbers never began to recover in these regions. Looking at abundance and biomass information from the survey there was evidence of strong recruitment among smaller individuals \< 1kg, but there has since not been any sizable population of fishes larger than 1kg outside of elasmobranchs. Work by [@pershing2015] suggests that part of the failure in recovery was due to an inability to account for temperature change in fisheries management. At this point in time the regional temperatures had just begun to reflect a regime shift, and could have been considered at that time an acute stressor. Since then the region has experienced nearly a decade of sustained above-average temperatures, and there are signs that the success seen in recruitment and survival of even the smaller size classes is declining.### Potential Drivers::::{.panel-tabset}#### Surface Temperature Anomalies ```{r}#| label: temperature-index#| fig.height: 6# Plot the timelinestemp_regimes %>%#filter(survey_area == "Full Region") %>% ggplot( aes(yr, sst_anom)) +#geom_col(aes(fill = regime), size = 0.75, alpha = 0.8) +geom_hline(yintercept =0, color =gmri_cols("orange"), size =1, lty =2) +geom_line(size =0.75, alpha =0.8) +scale_x_continuous(expand =expansion(add =c(0.25,0.25))) +scale_y_continuous(labels =number_format(suffix =" \u00b0C")) +facet_wrap(~survey_area, ncol =1) +guides(fill ="none") +theme(legend.title =element_blank(),legend.position ="bottom",legend.background =element_rect(fill ="transparent"),legend.key =element_rect(fill ="transparent", color ="transparent")) +labs(title ="NW-Atlantic Regional SST Anomalies",x ="", y ="Surface Temperature Anomaly",caption ="Anomalies calculated using 1982-2011 reference period.") +guides(label ="none",color =guide_legend(keyheight =unit(0.5, "cm")))```#### GARFO Landings```{r}#| label: landings-index#### Landings ####landings_summ %>%ggplot(aes(year, total_wt_z, group = survey_area)) +geom_hline(yintercept =0, color =gmri_cols("orange"), size =1, lty =2) +geom_line() +facet_wrap(~survey_area, ncol =1) +#scale_y_continuous(labels = scales::comma_format()) +scale_x_continuous(breaks =seq(1950, 2030, by =10)) +labs(x ="Year", y ="Finfish Landings Index (z)",title ="GARFO Finfish Landings")```#### Climate Modes```{r}#| label: climate-modes#### Climate Drivers ####ggplot(clim_drivers, aes(Time, Value, group = Var)) +geom_hline(yintercept =0, color =gmri_cols("orange"), size =1, lty =2) +geom_line() +facet_wrap(~Var, ncol =1) +scale_y_continuous(labels = scales::comma_format()) +scale_x_date(limits =as.Date(c("1965-01-01", "2020-12-31")),breaks = scales::breaks_pretty(n =6)) +labs(x ="Year", y ="Total Finfish Landings (lb.)",title ="Environmental Drivers")```#### CPR PCA```{r}#| label: cpr-modes#### CPR Community ##### CPR Plotcpr_indices %>%ggplot(aes(year, pca_loading, group = pca_mode)) +geom_hline(yintercept =0, color =gmri_cols("orange"), size =1, lty =2) +geom_line() +facet_wrap(~pca_mode, ncol =1) +scale_x_continuous(breaks =seq(1950, 2030, by =10)) +labs(x ="Year", y ="PCA Loading", title ="Gulf of Maine CPR PCA")```:::### All Correlations Matrix```{r}#| label: drivers-and-spectra#| eval: true#| fig-height: 7#| fig-width: 7# table to join for swapping shorthand for long-hand namesarea_df <-data.frame(area =c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "All"),survey_area =c("GoM", "GB", "SNE", "MAB", "Full Region"))# Common form for all:# year, area, var, value# Put climate drivers on annual scheduleclim_idx <- clim_drivers %>%group_by(year, area = EPU, var = Var) %>%summarise(value =mean(Value, na.rm = T),.groups ="drop")# landingslandings_idx <- landings_summ %>%group_by(year, area = survey_area) %>%summarise(value =mean(total_wt_z, na.rm = T),.groups ="drop") %>%mutate(var ="landings")# Zooplanktoncpr_idx <- cpr_indices %>%mutate(area ="Gulf of Maine",var =str_c("cpr_pca_", str_sub(pca_mode, -1, -1))) %>%group_by(year, area, var) %>%summarise(value =mean(pca_loading),.groups ="drop")# Temperaaturesst_idx <- temp_regimes %>%left_join(area_df) %>%select(year = yr, area, value = sst_anom) %>%mutate(var ="sst_anom")# Size Spectrum Slopes & Interceptsss_idx <- region_indices %>%select(year = Year, survey_area, spectra_slope = l10_slope_strat, spectra_int = l10_int_strat) %>%pivot_longer(cols =c(spectra_slope, spectra_int), names_to ="var", values_to ="value") %>%mutate(year =as.numeric(year)) %>%left_join(area_df) %>%select(-survey_area)# All driversall_drivers <-bind_rows(list( clim_idx, landings_idx, cpr_idx, sst_idx, ss_idx ))# Do some reshaping:driver_matrix <- all_drivers %>%mutate(id =str_replace_all(str_c(area, "_", var), " ", "_")) %>%select(-c(area, var)) %>%pivot_wider(names_from = id, values_from = value) %>%drop_na() %>%column_to_rownames(var ="year")# Correlation Plotrho <-cor(driver_matrix)corrplot(rho,method="ellipse",type="lower",order ="hclust",tl.col ="black",tl.srt =0,tl.cex =0.6,tl.offset =0.7,cl.cex =0.8,cl.offset =0.9,cl.ratio =0.1,title ="All Year Correlation Matrix")# GGplot correlation Matrix datagg_rho <-as.data.frame(rho) %>%rownames_to_column(var ="xvar") %>%pivot_longer(cols =-xvar, names_to ="yvar", values_to ="r2")``````{r}#| label: ggcorr-all#| fig-height: 6#| fig-width: 6#| eval: falsegg_rho %>%mutate(xvar =str_replace(xvar, "_", " "),yvar =str_replace(yvar, "_", " ")) %>%ggplot(aes(yvar, xvar, fill = r2)) +geom_tile() +scale_fill_distiller(palette ="RdBu", limits =c(-1, 1), breaks =c(-1,0,1)) +theme(axis.text.x =element_text(angle =45, hjust =1, vjust =1),strip.text.y =element_text(angle =0),legend.position ="bottom") +labs(x ="With these Other Indices (y)", y ="Correlation of These Indices (x)",fill ="Correlation Coefficient",title ="Correlation Matrix of All Drivers & Response")``````{r}#| label: ggcorr-drivers#| fig-height: 6#| fig-width: 8# Tidy it up, enhance the labelsgg_rho_tidy <- gg_rho %>%mutate(response_var =ifelse(str_detect(xvar, "slope|int"), T, F),driver_var =ifelse(str_detect(yvar, "osc|pca|landings|index|anom"), T, F),area =case_when(str_detect(xvar, "Georges") ~"Georges Bank",str_detect(xvar, "Gulf") ~"Gulf of Maine",str_detect(xvar, "Southern") ~"Southern New England",str_detect(xvar, "Mid-Atlantic") ~"Mid-Atlantic Bight"),area =factor(area, levels =c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")),driver_type =ifelse(str_detect(yvar, "landings"), "Fishing", "Environmental"),x_short =case_when(str_detect(xvar, "slope") ~"Spectra Slope",str_detect(xvar, "int") ~"Spectra Intercept"))# Do some filtering and string formatting:only_driver_response <- gg_rho_tidy %>%filter( response_var, driver_var) %>%mutate(yvar =str_replace_all(yvar, "_", " "),yvar =factor(yvar, levels =c("Gulf of Maine cpr pca 1","Gulf of Maine cpr pca 2","All gulf stream index","All north atlantic oscillation","All sst anom","Gulf of Maine sst anom","Georges Bank sst anom", "Southern New England sst anom","Mid-Atlantic Bight sst anom","Gulf of Maine landings","Georges Bank landings","Southern New England landings","Mid-Atlantic Bight landings" ))) # Plot the convoluted thing:only_driver_response %>%ggplot(aes(yvar, x_short, fill = r2)) +geom_tile() +facet_grid(area~driver_type, scales ="free_x", space ="free") +scale_fill_distiller(palette ="RdBu", limits =c(-1, 1), breaks =c(-1,0,1)) +theme(axis.text.x =element_text(angle =45, hjust =1, vjust =1),strip.text.y =element_text(angle =0),legend.position ="bottom",legend.title =element_text(vjust =0.75)) +labs(x ="Driver", y ="Bodymass Spectra Coefficient",fill ="Correlation Coefficient",title ="Correlation Between Size Spectra Response and Drivers",subtitle ="Correlation for Concurrent Years")```# References
